Appendix 3
Katie Florko
30/11/2022
1 Overview
This document contains the code from “A comparison of statistical models used to characterize species-habitat associations with movement data”. This code includes the case study analysis on ringed seal movement in eastern Hudson Bay, Canada.
2 Prepare
2.1 Load packages
library(tidyverse)
library(ggplot2)
library(lubridate)
library(rnaturalearth)
library(rgdal)
library(terra)
library(amt)
library(momentuHMM)
library(viridis)
library(sf)
library(here)2.2 Fish data
Next we load in the fish data. This is a subset of the data from Florko et al. (2021).
# load fish data
fish <- read.csv(here("data/preydiv.csv"))
head(fish)## lon lat preydiv
## 1 168329.9 -116109 0.513
## 2 168329.9 -166109 0.503
## 3 168329.9 -216109 0.512
## 4 168329.9 -266109 0.533
## 5 168329.9 -316109 0.557
## 6 168329.9 -366109 0.575Visualize the fish data.
#prepare world data
## List of azimuthal equal area - HudsonBay
natearth <- ne_countries(scale = "medium",returnclass = "sp")
natearth <- natearth[which(natearth$continent!="Antarctica"),]
nat_trans <- spTransform(natearth,CRS("+proj=laea +lat_0=60 +lon_0=-85 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"))## Warning in spTransform(xSP, CRSobj, ...): NULL source CRS comment, falling back
## to PROJ string## Warning in wkt(obj): CRS object has no comment# plot fish map
fishmap <- ggplot() +
geom_tile(data = fish, aes(x = lon,y = lat,fill = preydiv)) +
scale_fill_viridis(option = "mako", limits = c(0.5, 0.75), name = "Prey\ndiversity") +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F)
fishmap
2.3 Seal data
Next we load in the seal movement data. This is a subset of the data from Florko et al. (2023).
seal <- read.csv(here("data/seal_track_m.csv"))
head(seal)## lon lat date id
## 1 357940.9 -368949.6 2012-10-29 14:00 116484_1
## 2 371594.3 -353092.9 2012-10-30 14:00 116484_1
## 3 388429.2 -361478.4 2012-10-31 14:00 116484_1
## 4 399326.2 -357288.1 2012-11-01 14:00 116484_1
## 5 411788.2 -349500.7 2012-11-02 14:00 116484_1
## 6 407708.1 -372564.6 2012-11-03 14:00 116484_1# ensure the data is in the correct format
seal <- seal %>%
mutate(id = as.character(id),
date = as.Date(date)) Visualize the seal data on top of the fish data.
# plot
sealfishmap <- fishmap +
geom_point(data=seal, aes(x=lon, y=lat), alpha = 0.6, color = "#FCEEAE") +
geom_path(data=seal, aes(x=lon, y=lat), color = "#FCEEAE")
sealfishmap
3 Fit models
3.1 RSF
Fit using amt (Signer et al. 2019)
# generate availability sample
set.seed(2023)
data_rsf <- seal %>%
make_track(lon, lat, date) %>%
random_points() # default is ten times as many available points as observed points
# plot used vs available locations
sealfishmap +
geom_point(data=data_rsf, aes(x=x_, y=y_, color = case_), alpha = 0.6) +
scale_color_manual(values = c("black", "#FCEEAE"),
label = c("Available", "Observed"), name = "Data type") 
# rasterize and extract prey diversity covariate
fish_raster <- terra::rast(fish, crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
data_rsf <- data_rsf %>%
extract_covariates(fish_raster)
# fit rsf
rsf1 <- data_rsf %>%
amt::fit_rsf(case_ ~ preydiv, model = TRUE)
# see parameter estimates
summary(rsf1)##
## Call:
## stats::glm(formula = formula, family = stats::binomial(link = "logit"),
## data = data, model = TRUE)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.5753 -0.4729 -0.4294 -0.3754 2.3363
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.118 1.400 -4.369 1.25e-05 ***
## preydiv 6.353 2.312 2.748 0.006 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 938.28 on 1539 degrees of freedom
## Residual deviance: 930.60 on 1538 degrees of freedom
## AIC: 934.6
##
## Number of Fisher Scoring iterations: 53.2 SSF
Fit using amt (Signer et al. 2019)
# generate availability sample
set.seed(2023)
data_ssf <- seal %>%
make_track(lon, lat, date) %>%
steps() %>%
random_steps() %>% # default is ten times as many available points as observed points
arrange(case_)
# plot used vs available locations
sealfishmap +
geom_point(data=data_ssf, aes(x=x2_, y=y2_, color = case_), alpha = 0.6) +
scale_color_manual(values = c("black", "#FCEEAE"),
label = c("Available", "Observed"), name = "Data type") 
# extract prey diversity covariate
data_ssf <- data_ssf %>%
extract_covariates(fish_raster, where = "end") #sample at end of step
# transform movement covariates
data_ssf <- data_ssf %>%
mutate(cos_ta = cos(ta_),
log_sl = log(sl_))
# fit ssfs
## ssf1: ssf without covariate affecting movement kernel
ssf1 <- data_ssf %>%
fit_clogit(case_ ~ preydiv + log_sl + cos_ta + strata(step_id_), model = TRUE)
## ssf2: ssf with covariate affecting movement kernel
ssf2 <- data_ssf %>%
fit_clogit(case_ ~ log_sl * preydiv + cos_ta + strata(step_id_), model = TRUE)
# see parameter estimates
summary(ssf1)## Call:
## coxph(formula = Surv(rep(1, 1518L), case_) ~ preydiv + log_sl +
## cos_ta + strata(step_id_), data = data, model = TRUE, method = "exact")
##
## n= 1516, number of events= 138
## (2 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## preydiv 0.957799 2.605953 6.772711 0.141 0.888
## log_sl -0.008404 0.991631 0.083328 -0.101 0.920
## cos_ta 0.031707 1.032215 0.130643 0.243 0.808
##
## exp(coef) exp(-coef) lower .95 upper .95
## preydiv 2.6060 0.3837 4.477e-06 1.517e+06
## log_sl 0.9916 1.0084 8.422e-01 1.168e+00
## cos_ta 1.0322 0.9688 7.990e-01 1.333e+00
##
## Concordance= 0.494 (se = 0.029 )
## Likelihood ratio test= 0.09 on 3 df, p=1
## Wald test = 0.09 on 3 df, p=1
## Score (logrank) test = 0.09 on 3 df, p=1summary(ssf2)## Call:
## coxph(formula = Surv(rep(1, 1518L), case_) ~ log_sl * preydiv +
## cos_ta + strata(step_id_), data = data, model = TRUE, method = "exact")
##
## n= 1516, number of events= 138
## (2 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## log_sl 1.663e+00 5.273e+00 1.286e+00 1.293 0.196
## preydiv 2.791e+01 1.318e+12 2.177e+01 1.282 0.200
## cos_ta 2.338e-02 1.024e+00 1.311e-01 0.178 0.858
## log_sl:preydiv -2.744e+00 6.433e-02 2.101e+00 -1.306 0.192
##
## exp(coef) exp(-coef) lower .95 upper .95
## log_sl 5.273e+00 1.897e-01 4.241e-01 6.556e+01
## preydiv 1.318e+12 7.589e-13 3.919e-07 4.430e+30
## cos_ta 1.024e+00 9.769e-01 7.917e-01 1.324e+00
## log_sl:preydiv 6.433e-02 1.555e+01 1.046e-03 3.955e+00
##
## Concordance= 0.555 (se = 0.029 )
## Likelihood ratio test= 1.81 on 4 df, p=0.8
## Wald test = 1.81 on 4 df, p=0.8
## Score (logrank) test = 1.81 on 4 df, p=0.83.3 HMM
Fit using momentuHMM (McClintock and Michelot, 2018)
# prepare dataset
data_hmm <- seal %>%
make_track(lon, lat, date) %>%
extract_covariates(fish_raster) %>%
mutate(ID = 1,
x = x_,
y = y_,
date = t_) %>%
dplyr::select(ID, x, y, date, preydiv) %>%
as.data.frame()
# convert into momentuHMM format
data_hmm <- momentuHMM::prepData(data_hmm, coordNames = c("x", "y"),
type = "UTM",
covNames = c("preydiv"))
# define initial parameters
nbStates <- 3 # number of states
stepDist <- "gamma" # step distribution
angleDist <- "vm" # turning angle distribution
mu0 <- c(5000, 18000, 38000) # mean step length for each state
sigma0 <- c(5000, 10000, 10000) # sd step length for each state
stepPar0 <- c(mu0, sigma0)
kappa0 <- c(0.35, 0.55, 0.5) # turning angle for each state
formula = ~ preydiv # identify covariates
# fit basic 3-state hmm
set.seed(2023)
hmm1 <- momentuHMM::fitHMM(data=data_hmm, nbStates=nbStates,
dist=list(step=stepDist,angle=angleDist),
Par0=list(step=stepPar0,angle=kappa0))
# retrieve parameters to refine model
Par0_hmm1 <- momentuHMM::getPar0(hmm1, formula=formula)
# fit a refined hmm with parameters from hhm1
set.seed(2023)
hmm2 <- momentuHMM::fitHMM(data=data_hmm, nbStates=3,
dist=list(step=stepDist,angle=angleDist),
Par0=Par0_hmm1$Par,
delta0 = Par0_hmm1$delta,
beta0 = Par0_hmm1$beta,
formula=formula)
# plot decoded states
data_hmm$state <- as.factor(momentuHMM::viterbi(hmm2))
hmmstate_plot <- ggplot() +
scale_fill_viridis(option = "mako", limits = c(0.5, 0.75), name = "Prey\ndiversity") +
geom_path(data=data_hmm, aes(x=x, y=y, color = state, group =ID)) +
geom_point(data=data_hmm, aes(x=x, y=y, color = state, shape = state), size=2, alpha = 0.8) +
scale_color_manual(values = c("#99DDB6", "#539D9C", "#312C66"),
labels = c("Localized\nsearch", "Area-restricted\nsearch", "Travelling"),
name = "HMM state") +
scale_shape_manual(values = c(15, 16, 17),
labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling"),
name = "HMM state") +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F)
hmmstate_plot
#---------------------PLOT HISTOGRAMS
## prep data in lat/lon and refit model (so step length is in kilometers)
data_hmm <- seal %>%
make_track(lon, lat, date) %>%
extract_covariates(fish_raster) %>%
mutate(ID = 1,
x = x_,
y = y_,
date = t_) %>%
dplyr::select(ID, x, y, date, preydiv) %>%
as.data.frame() %>%
st_as_sf(coords = c("x", "y")) %>%
st_set_crs("+proj=laea +lat_0=60 +lon_0=-85 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0") %>%
st_transform("+proj=longlat +datum=WGS84") %>%
mutate(long = unlist(map(geometry,1)),
lat = unlist(map(geometry,2))) %>%
st_drop_geometry()
# do momentuhmm which calculates step length in KM
data_hmm <- momentuHMM::prepData(data_hmm, coordNames = c("long", "lat"), type = "LL", covNames = c("preydiv"))
# define parameters
nbStates <- 3
stepDist <- "gamma" # step distribution
angleDist <- "vm" # turning angle distribution
mu0 <- c(5, 12, 38)
sigma0 <- c(3, 5, 8)
stepPar0 <- c(mu0, sigma0)
kappa0 <- c(0.35, 0.55, 0.5)
# fit HMM (with step in km)
set.seed(2023)
hmm1_km <- momentuHMM::fitHMM(data=data_hmm, nbStates=nbStates,
dist=list(step=stepDist,angle=angleDist),
Par0=list(step=stepPar0,angle=kappa0))
# get parameters
Par0_hmm1_km <- momentuHMM::getPar0(hmm1_km, formula=formula)
# fit HMM with parameters from hmm1_km
set.seed(2023)
hmm2_km <- momentuHMM::fitHMM(data=data_hmm, nbStates=3,
dist=list(step=stepDist,angle=angleDist),
Par0=Par0_hmm1_km$Par,
delta0 = Par0_hmm1_km$delta,
beta0 = Par0_hmm1_km$beta,
formula=formula)
# add the state estimate from the HMM
data_hmm$state <- as.factor(momentuHMM::viterbi(hmm2_km))
##-----STEP DIST histogram
#calculate frequencies of states
v <- momentuHMM::viterbi(hmm2_km)
stateFreq <- table(v) / length(v)
#plot colours
colours.states <- c("#99DDB6", "#539D9C", "#312C66")
#generate sequence for x axis of density functions
x <- seq(0, 50, length=1000)
#get converged mean and sd for each state
meanARS <- hmm2_km$mle$step[1,1]
sdARS <- hmm2_km$mle$step[2,1]
meanCR <- hmm2_km$mle$step[1,2]
sdCR <- hmm2_km$mle$step[2,2]
meanTR <- hmm2_km$mle$step[1,3]
sdTR <- hmm2_km$mle$step[2,3]
#calculate shape and scale of the gamma distributions from mean and sd
sh <- function(mean, sd) { return(mean^2 / sd^2)}
sc <- function(mean, sd) { return(sd^2 / mean)}
#get density functions of the distributions
y_ARS <- dgamma(x, shape=sh(meanARS,sdARS), scale=sc(meanARS,sdARS)) * stateFreq[[1]]
y_CR <- dgamma(x, shape=sh(meanCR,sdCR), scale=sc(meanCR,sdCR)) * stateFreq[[2]]
y_TR <- dgamma(x, shape=sh(meanTR,sdTR), scale=sc(meanTR,sdTR)) * stateFreq[[3]]
#combine densities in a single dataframe for more convenient plotting
df.y_ARS <- data.frame(dens=y_ARS, State="Foraging", x=x)
df.y_CR <- data.frame(dens=y_CR, State="ARS", x=x)
df.y_TR <- data.frame(dens=y_TR, State="Travelling", x=x)
statedis <- rbind(df.y_ARS, df.y_CR, df.y_TR)
#plot distributions
hmmstepdist_plot <- ggplot() +
geom_line(data=statedis,aes(x=x,y=dens,colour=State,linetype=State), size=1.2) +
scale_colour_manual(values=c(colours.states,"#000000"),
breaks = c('Foraging', 'ARS', 'Travelling'),
labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling")) +
scale_linetype_manual(values=c("solid","solid", "solid"),
breaks = c('Foraging', 'ARS', 'Travelling'),
labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling")) +
ylab("Density") +
xlab("Step length (kms)") +
theme_minimal()## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## Warning: Please use `linewidth` instead.hmmstepdist_plot
##-----ANGLE DIST
# #generate sequence for x axis of density functions
x <- seq(-pi, pi,length=1000)
#get converged mean and concentration for each state
meanARS <- hmm2_km$mle$angle[1,1]
sdARS <- hmm2_km$mle$angle[2,1]
meanCR <- hmm2_km$mle$angle[1,2]
sdCR <- hmm2_km$mle$angle[2,2]
meanTR <- hmm2_km$mle$angle[1,3]
sdTR <- hmm2_km$mle$angle[2,3]
#get density functions of the distributions
y_ARS <- CircStats::dvm(x, mu=meanARS, kappa=sdARS) * stateFreq[[1]]
y_CR <- CircStats::dvm(x, mu=meanCR, kappa=sdCR) * stateFreq[[2]]
y_TR <- CircStats::dvm(x, mu=meanTR, kappa=sdTR) * stateFreq[[3]]
#combine densities in a single dataframe for more convenient plotting
df.y_ARS <- data.frame(dens=y_ARS, State="Foraging",x=x)
df.y_CR <- data.frame(dens=y_CR, State="ARS",x=x)
df.y_TR <- data.frame(dens=y_TR, State="Travelling",x=x)
cmb <- rbind(df.y_TR, df.y_CR, df.y_ARS)
#plot distributions
hmmangledist_plot <- ggplot() +
geom_line(data=cmb,aes(x=x,y=dens,colour=State), size = 1.2) +
scale_colour_manual(values=c(colours.states), breaks = c('Foraging', 'ARS', 'Travelling'), labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travelling")) +
scale_x_continuous(limits=c(-pi,pi))+
ylab("Density") +
xlab("Turn angle (radians) ") +
theme_minimal()
hmmangledist_plot
4 Plot prediction maps
Prep the data
# prep the fish data
newfish <- fish_raster %>%
terra::as.data.frame(xy = TRUE) %>%
filter(x > 100000 & x < 600000 & y > -550000 & y < 0)4.1 RSF
modcoef <- summary(rsf1)$coef
x <- exp(modcoef[1] + (modcoef[2] * newfish$preydiv))
x <- x / (1 + x)
# range fn
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
newfish$rsf_prediction <-range01(x)
# plot
map_RSF <- ggplot() +
geom_tile(data = newfish, aes(x = x,y = y, fill = rsf_prediction)) +
scale_fill_viridis(option = "mako", name = "RSF prediction", limits = c(0, 0.3)) +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) ## Regions defined for each Polygonsmap_RSF
4.2 SSF
# generate availability sample
set.seed(2023)
data_ssf <- seal %>%
mutate(date = as.POSIXct(date)) %>%
make_track(lon, lat, date) %>%
steps() %>%
random_steps() %>%
arrange(case_) %>%
amt::extract_covariates(fish_raster, where = "both") %>%
na.omit()
### -------------- SSF1
# fit SSF1 model
m1 <- data_ssf |>
fit_clogit(case_ ~ preydiv_end + cos(ta_) + log(sl_) +
# add terms for home range
x2_ + y2_ + I(x2_^2 + y2_^2) +
strata(step_id_))
# Start with simulations
# First generate a redistribution kernel
set.seed(2023)
start <- make_start((seal %>%
mutate(date = as.POSIXct(date)) %>%
make_track(lon, lat, date))
# random sample along the track as the start
[sample(1:nrow(seal), 1), ],
dt_ = hours(2))
# generate redistribution kernel
k1 <- redistribution_kernel(m1, map = fish_raster, start = start,
stochastic = TRUE,
tolerance.outside = 0.1,
n.control = 1e3)
# Now simulate a path of length 30 for 50 times. This will lead a transient UD. For a steady state UD, either increase `n` or choose many more random starting points.
n_steps = 1e3
n_steps1 = n_steps + 1
set.seed(2023)
p1 <- amt::simulate_path(x = k1, n = n_steps, start = start, verbose = TRUE)
# plot simulated track
ssf_track_1 <- fishmap +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
geom_point(data = p1, aes(x = x_,y = y_), alpha = 0.61) +
geom_path(data = p1, aes(x = x_,y = y_)) ## Regions defined for each Polygons# estimate SSUD
uds_ssf1 <- tibble(rep = 1:n_steps1,
x_ = p1$x_, y_ = p1$y_,
t_ = p1$t_, dt = p1$dt) |>
filter(!is.na(x_)) |>
make_track(x_, y_) |>
hr_kde(trast = fish_raster, which_min = "global") %>%
hr_ud() %>%
terra::as.data.frame(xy = TRUE)
# plot SSUD
map_SSF1 <- ggplot() +
geom_tile(data = uds_ssf1, aes(x = x,y = y, fill = preydiv)) +
scale_fill_viridis(option = "mako", name = "SSF1 prediction") +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) ## Regions defined for each Polygonsmap_SSF1
### -------------- SSF2
# fit SSF2 model
m2 <- data_ssf |>
fit_clogit(case_ ~ preydiv_end + preydiv_end:cos(ta_)+
preydiv_end:log(sl_) +
# add terms for home range
x2_ + y2_ + I(x2_^2 + y2_^2) +
strata(step_id_))
# Start with simulations
# First generate a redistribution kernel
set.seed(2023)
start <- make_start((seal %>%
mutate(date = as.POSIXct(date)) %>%
make_track(lon, lat, date))
# first location of the track as the start
[sample(1:nrow(seal), 1), ],
dt_ = hours(2))
# generate redistribution kernel
k2 <- redistribution_kernel(m2, map = fish_raster, start = start,
stochastic = TRUE,
tolerance.outside = 0.3,
n.control = 1e3)
# Now simulate a path of length 30 for 50 times. This will lead a transient UD. For a steady state UD, either increase `n` or choose many more random starting points.
n_steps = 1e3
n_steps1 = n_steps + 1
set.seed(2023)
p2 <- amt::simulate_path(x = k2, n = n_steps, start = start)
# plot simulated track
ssf_track_2 <- fishmap +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
geom_point(data = p2, aes(x = x_,y = y_), alpha = 0.61) +
geom_path(data = p2, aes(x = x_,y = y_)) ## Regions defined for each Polygons# estimate SSUD
uds_ssf2 <- tibble(rep = 1:n_steps1,
x_ = p1$x_, y_ = p1$y_,
t_ = p1$t_, dt = p1$dt) |>
filter(!is.na(x_)) |>
make_track(x_, y_) |>
hr_kde(trast = fish_raster, which_min = "local") %>%
hr_ud() %>%
terra::as.data.frame(xy = TRUE)
# plot SSUD
map_SSF2 <- ggplot() +
geom_tile(data = uds_ssf2, aes(x = x,y = y, fill = preydiv)) +
scale_fill_viridis(option = "mako", name = "SSF2 prediction") +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), fill = "grey80", color = "white") +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) ## Regions defined for each Polygonsmap_SSF2
library(cowplot)##
## Attaching package: 'cowplot'
##
## The following object is masked from 'package:lubridate':
##
## stampplot_grid(map_SSF1, map_SSF2, ssf_track_1, ssf_track_2, nrow =2)
4.3 HMM
x <- as.data.frame(momentuHMM::stationary(hmm2, data.frame(preydiv = newfish$preydiv)))
newfish$hmm_state1 <- x$state.1
newfish$hmm_state2 <- x$state.2
newfish$hmm_state3 <- x$state.3
newfish_long <- newfish %>%
tidyr::pivot_longer(cols = hmm_state1:hmm_state3,
names_to = "model", values_to = "prediction") %>%
mutate(dplyr::across(model, factor, levels=
c("hmm_state1", "hmm_state2", "hmm_state3")))
map_HMM <- ggplot() +
geom_tile(data = newfish_long, aes(x = x, y = y, fill = prediction)) +
scale_fill_viridis(option = "mako", limits = c(0,1)) +
labs(fill = 'HMM predicted\nprobability') +
geom_polygon(data = nat_trans, aes(x=long,y=lat,group=group), colour = "white", fill = "grey90", size = 0.3) +
coord_cartesian(xlim = c(150000,550000), ylim = c(-500000,-50000), expand = F) +
facet_wrap(~model, labeller = as_labeller(c('hmm_state1' = "Localized search",
'hmm_state2' = "Area-restricted search",
'hmm_state3' = "Travel")))## Regions defined for each Polygonsmap_HMM
5 Plot prey div vs probability of selection
5.1 RSF
base <- newfish %>%
# (note, it doesn't matter what values you pick for sl_ and ta_)
mutate(log_sl = log(45),
cos_ta = cos(1))
x1 <- base
x2 <- base %>%
mutate(preydiv = mean(base$preydiv))
# So how do we iterate this? We could use a for loop, we could use lapply, etc.
# I will use lapply here, with a custom function inside.
log_rss_list <- lapply(1:nrow(x1), function(i) {
# Report status
cat(i, "of", nrow(x1), "\n")
# Calculate log-RSS for that row
xx <- log_rss(rsf1, x1[i,], x2[i,], ci = "se")
# Return the element $df
return(xx$df)
})## 1 of 72
## 2 of 72
## 3 of 72
## 4 of 72
## 5 of 72
## 6 of 72
## 7 of 72
## 8 of 72
## 9 of 72
## 10 of 72
## 11 of 72
## 12 of 72
## 13 of 72
## 14 of 72
## 15 of 72
## 16 of 72
## 17 of 72
## 18 of 72
## 19 of 72
## 20 of 72
## 21 of 72
## 22 of 72
## 23 of 72
## 24 of 72
## 25 of 72
## 26 of 72
## 27 of 72
## 28 of 72
## 29 of 72
## 30 of 72
## 31 of 72
## 32 of 72
## 33 of 72
## 34 of 72
## 35 of 72
## 36 of 72
## 37 of 72
## 38 of 72
## 39 of 72
## 40 of 72
## 41 of 72
## 42 of 72
## 43 of 72
## 44 of 72
## 45 of 72
## 46 of 72
## 47 of 72
## 48 of 72
## 49 of 72
## 50 of 72
## 51 of 72
## 52 of 72
## 53 of 72
## 54 of 72
## 55 of 72
## 56 of 72
## 57 of 72
## 58 of 72
## 59 of 72
## 60 of 72
## 61 of 72
## 62 of 72
## 63 of 72
## 64 of 72
## 65 of 72
## 66 of 72
## 67 of 72
## 68 of 72
## 69 of 72
## 70 of 72
## 71 of 72
## 72 of 72# We now have a list with an element for each row. We can combine those rows
# into a single data.frame with 'dplyr::bind_rows()'
res <- dplyr::bind_rows(log_rss_list)
# PLOT
rsf_line <- ggplot(res, aes(x = preydiv_x1, y = (log_rss))) +
geom_line(size = 1, color = "#F8D59F",linetype = 2) +
geom_ribbon(aes(ymin=lwr, ymax=upr, x=preydiv_x1), alpha = 0.4, fill = "#F8D59F") +
xlab("Prey diversity") +
ylab("log-RSS") +
theme_minimal()5.2 SSF
## Prediction for SSF1
x1 <- base
x2 <- base %>%
mutate(preydiv = mean(base$preydiv))
# So how do we iterate this? We could use a for loop, we could use lapply, etc.
# I will use lapply here, with a custom function inside.
log_rss_list <- lapply(1:nrow(x1), function(i) {
# Calculate log-RSS for that row
xx <- log_rss(ssf1, x1[i,], x2[i,], ci = "se")
# Return the element $df
return(xx$df)
})
# We now have a list with an element for each row. We can combine those rows
# into a single data.frame with 'dplyr::bind_rows()'
res <- dplyr::bind_rows(log_rss_list)
# PLOT
ssf_line_1 <- ggplot() +
geom_line(data = res, aes(x = preydiv_x1, y = log_rss), size = 1, linetype = 3) +
geom_ribbon(data = res, aes(ymin=lwr, ymax=upr, x=preydiv_x1), alpha = 0.3) +
xlab("Prey diversity") +
ylab("log-RSS") +
theme_minimal()
## Prediction for SSF2
nums = seal %>%
make_track(lon, lat, date) %>%
steps %>%
summarize(quants = quantile(sl_, c(0.25, 0.5, 0.75))) %>%
pull()
x1 <- base
results <- lapply(nums, function(i) {
x1$log_sl <- i
x2 <- base %>%
mutate(preydiv = mean(base$preydiv))
# So how do we iterate this? We could use a for loop, we could use lapply, etc.
# I will use lapply here, with a custom function inside.
log_rss_list <- lapply(1:nrow(x1), function(i) {
# Calculate log-RSS for that row
xx <- log_rss(ssf2, x1[i,], x2[i,], ci = "se")
# Return the element $df
return(xx$df)
})
# We now have a list with an element for each row. We can combine those rows
# into a single data.frame with 'dplyr::bind_rows()'
res <- dplyr::bind_rows(log_rss_list)
} )
results <- dplyr::bind_rows(results) %>%
mutate(log_sl_x1 = as.factor(round(log_sl_x1, 1)))
# PLOT
ggplot(results, aes(x = preydiv_x1, y = (log_rss))) +
geom_line(size = 1, aes(color = log_sl_x1, group = log_sl_x1)) +
geom_ribbon(aes(ymin=lwr, ymax=upr, x=preydiv_x1, fill = log_sl_x1, group = log_sl_x1), alpha = 0.3) +
xlab("Prey diversity") +
ylab("log-RSS") +
theme_minimal()
# Plot
ssf_line_1 +
geom_line(data = results, size = 1, aes(x = preydiv_x1, y = (log_rss),
color = log_sl_x1, group = log_sl_x1)) +
geom_ribbon(data = results, aes(ymin=lwr, ymax=upr, x=preydiv_x1,
fill = log_sl_x1, group = log_sl_x1), alpha = 0.3)
5.3 HMM
# grab stationary probabilities
ps <- momentuHMM::plotStationary(hmm2, plotCI= TRUE, return = TRUE)
state1 <- ps$preydiv$'state 1' %>% mutate(state = 1)
state2 <- ps$preydiv$'state 2' %>% mutate(state = 2)
state3 <- ps$preydiv$'state 3' %>% mutate(state = 3)
# bind to one data frame
pdat <- rbind(state1, state2, state3) %>%
mutate(state = as.character(state)) # plot
hmmprobs_plot <- ggplot() +
geom_line(data = pdat, aes(x = cov, y = est, color = state)) +
geom_ribbon(data = pdat, aes(x=cov, y=est, ymax=est+se, ymin=est-se, fill = state),
alpha = 0.4, show.legend = TRUE) +
ylab("Stationary state probabilties") +
xlab("Prey diversity") +
scale_color_manual(values = c("#99DDB6", "#539D9C", "#312C66"), name = "HMM state",
labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travel")) +
scale_fill_manual(values = c("#99DDB6", "#539D9C", "#312C66"), name = "HMM state",
labels=c("Localized\nsearch", "Area-restricted\nsearch", "Travel")) +
theme_minimal()
hmmprobs_plot6 References
Florko, K.R.N., Tai, T.C., Cheung, W.W.L., Sumaila, U.R., Ferguson, S.H., Yurkowski, D.J., Auger-Méthé, M. 2021. Predicting how climate change threatens the prey base of Arctic marine predators. Ecology Letters, 24: 2563-2575. https://doi.org/10.1111/ele.13866
Florko, K.R.N., Shuert, C.R., Cheung, W.W.L., Ferguson, S.H., Jonsen, I.D., Rosen, D.A.S., Sumaila, U.R., Tai, T.C., Yurkowski, D.J., Auger-Méthé, M. 2023. Linking movement and dive data to prey distribution models: new insights in foraging behavior and potential pitfalls of movement analyses. Movement Ecology, 11:17 https://doi.org/10.1186/s40462-023-00377-2
McClintock, B.T., Michelot, T. 2018. momentuHMM: R package for generalized hidden markov models of animal movement. Methods Ecol. Evolut. 9, 1518-1530. https://doi.org/10.1111/2041-210X.12995
Signer, J., Fieberg, J., & Avgar, T. (2017). Estimating utilization distributions from fitted step‐selection functions. Ecosphere, 8(4), e01771. https://doi.org/10.1002/ecs2.1771
Signer, J., J. Fieberg, and T. Avgar. 2019. Animal movement tools (amt): R package for managing tracking data and conducting habitat selection analyses. Ecology and Evolution 9:880–890. https://doi.org/10.1002/ece3.4823